#ifndef RECTCELLS_HPP
#define RECTCELLS_HPP

#include "structCellsBase.hpp"
#include "eigenConcat.hpp"

namespace hms
{

namespace rectCells_internal
{

using AreasExpression =
	Eigen::Reshaped<const Eigen::ArrayWrapper<const Eigen::Product<
		Eigen::MatrixWrapper<const ArrayXs >,
		Eigen::Transpose<const Eigen::MatrixWrapper<const ArrayXs>>,
		0>>,
	-1, 1, 0>;

using CircumfsExpression  = const Eigen::CwiseBinaryOp<
	Eigen::internal::scalar_product_op<scalar, scalar>,
	const Eigen::CwiseNullaryOp<
		Eigen::internal::scalar_constant_op<scalar>, const ArrayXs
	>,
	const Eigen::Reshaped<
		const Eigen::CwiseBinaryOp<
			Eigen::internal::scalar_sum_op<scalar, scalar>,
			const Eigen::Replicate<ArrayXs, -1, -1>,
			const Eigen::Replicate<Eigen::Transpose<const ArrayXs>,-1, -1>
		>,
		-1, 1, 0
	> >;

using ArraySegment = const Eigen::Block<const ArrayXs, -1, 1, false>;
/* sum of neighbouring coordinates divided by 2 (=middle between them) */
// using CentroidCoords = 
// 	Eigen::MatrixWrapper<const Eigen::CwiseBinaryOp<
// 		Eigen::internal::scalar_product_op<scalar, scalar>,
// 		const Eigen::CwiseNullaryOp<
// 			Eigen::internal::scalar_constant_op<scalar>,
// 			const ArrayXs
// 		>,
// 		const Eigen::CwiseBinaryOp<
// 			Eigen::internal::scalar_sum_op<scalar, scalar>,
// 			ArraySegment,
// 			ArraySegment
// 		>
// 	> >;
using CentroidCoords = const Eigen::MatrixWrapper<const ArrayXs>;

using CentroidsExpression = EigenCustom::ConcatReturnType<
	EigenCustom::Direction::vertical,
	Eigen::Replicate<Eigen::Transpose<CentroidCoords>, -1, -1>,
	Eigen::Transpose<Eigen::Reshaped<
		const Eigen::Replicate<Eigen::Transpose<CentroidCoords>,-1, -1>,
		-1, 1, 0
	> > >;

} // namespace rectCells_internal


template<>
class CellsTraits<RectCells>
{
public:
	using MeshType = RectMesh;
	using AreasReturnType =
		#ifdef RECTMESH_PRECALC_AREAS
			const ArrayXs&;
		#else
			rectCells_internal::AreasExpression;
		#endif

	using CircumfsReturnType = 
		#ifdef RECTMESH_PRECALC_CIRCUMFS
			const ArrayXs&;
		#else
			rectCells_internal::CircumfsExpression;
		#endif

	using CentroidsReturnType =
		#ifdef RECTMESH_PRECALC_CENTROIDS
			const Matrix2Xs&;
		#else
			rectCells_internal::CentroidsExpression;
		#endif
};


class RectCells : public StructCellsBase<RectCells>
{
public:
	using Traits = CellsTraits<RectCells>;
	using AreasReturnType = typename Traits::AreasReturnType;
	using CircumfsReturnType = typename Traits::CircumfsReturnType;
	using CentroidsReturnType = typename Traits::CentroidsReturnType;

private:
	ArrayXs m_xCoords;
	ArrayXs m_yCoords;

	ArrayXs m_dx;
	ArrayXs m_dy;

#ifdef RECTMESH_PRECALC_AREAS
	ArrayXs m_areas;
#endif

#ifdef RECTMESH_PRECALC_CIRCUMFS
	ArrayXs m_circumfs;
#endif

#ifdef RECTMESH_PRECALC_CENTROIDS
	Matrix2Xs m_centroids;
#endif

public:
	RectCells() = default;
	RectCells( const RectMesh& );

	auto xCoords() const -> const ArrayXs&;
	auto yCoords() const -> const ArrayXs&;

	auto dx() const -> const ArrayXs&;
	auto dy() const -> const ArrayXs&;
private:
	auto areasImp    () const -> AreasReturnType;
	auto circumfsImp () const -> CircumfsReturnType;
	auto centroidsImp() const -> CentroidsReturnType;

	void storeData();

	friend CellsBase;
	friend RectMesh;
};

} // namespace hms

#endif